Including a spatial predictive process in band recovery models improves inference for Lincoln estimates of animal abundance

Abstract Abundance estimation is a critical component of conservation planning, particularly for exploited species where managers set regulations to restrict harvest based on current population size. An increasingly common approach for abundance estimation is through integrated population modeling (IPM), which uses multiple data sources in a joint likelihood to estimate abundance and additional demographic parameters. Lincoln estimators are one commonly used IPM component for harvested species, which combine information on the rate and total number of individuals harvested within an integrated band‐recovery framework to estimate abundance at large scales. A major assumption of the Lincoln estimator is that banding and recoveries are representative of the whole population, which may be violated if major sources of spatial heterogeneity in survival or harvest rates are not incorporated into the model. We developed an approach to account for spatial variation in harvest rates using a spatial predictive process, which we incorporated into a Lincoln estimator IPM. We simulated data under different configurations of sample sizes, harvest rates, and sources of spatial heterogeneity in harvest rate to assess potential model bias in parameter estimates. We then applied the model to data collected from a field study of wild turkeys (Meleagris gallapavo) to estimate local and statewide abundance in Maine, USA. We found that the band recovery model that incorporated a spatial predictive process consistently provided estimates of adult and juvenile abundance with low bias across a variety of spatial configurations of harvest rate and sampling intensities. When applied to data collected on wild turkeys, a model that did not incorporate spatial heterogeneity underestimated the harvest rate in some subregions. Consistent with simulation results, this led to overestimation of both local and statewide abundance. Our work demonstrates that a spatial predictive process is a viable mechanism to account for spatial variation in harvest rates and limit bias in abundance estimates. This approach could be extended to large‐scale band recovery data sets and has applicability for the estimation of population parameters in other ecological models as well.

Lincoln estimators are one commonly used IPM component for harvested species, which combine information on the rate and total number of individuals harvested within an integrated band-recovery framework to estimate abundance at large scales.
A major assumption of the Lincoln estimator is that banding and recoveries are representative of the whole population, which may be violated if major sources of spatial heterogeneity in survival or harvest rates are not incorporated into the model. We developed an approach to account for spatial variation in harvest rates using a spatial predictive process, which we incorporated into a Lincoln estimator IPM. We simulated data under different configurations of sample sizes, harvest rates, and sources of spatial heterogeneity in harvest rate to assess potential model bias in parameter estimates. We then applied the model to data collected from a field study of wild turkeys (Meleagris gallapavo) to estimate local and statewide abundance in Maine, USA.
We found that the band recovery model that incorporated a spatial predictive process consistently provided estimates of adult and juvenile abundance with low bias across a variety of spatial configurations of harvest rate and sampling intensities. When applied to data collected on wild turkeys, a model that did not incorporate spatial heterogeneity underestimated the harvest rate in some subregions. Consistent with simulation results, this led to overestimation of both local and statewide abundance.
Our work demonstrates that a spatial predictive process is a viable mechanism to account for spatial variation in harvest rates and limit bias in abundance estimates. This approach could be extended to large-scale band recovery data sets and has applicability for the estimation of population parameters in other ecological models as well.

| INTRODUC TI ON
Abundance estimation is a critical component of successful conservation planning (Thogmartin et al., 2006), particularly for exploited species where managers set quotas or other regulations to restrict harvest based on current population size (Nichols et al., 2007;Runge et al., 2009). If abundance is overestimated, liberal regulations may lead to a larger portion of the population being removed than intended, which can negatively impact long-term stability (Johnson et al., 2012;Weinbaum et al., 2013). Alternatively, if population size is underestimated, harvest regulations may be set more restrictively than necessary, leading to underutilization of the resource, and reducing the opportunity for consumers. In either instance, there are benefits to identifying and implementing tools that estimate abundance as accurately as possible.
An increasingly common approach for abundance estimation is integrated population modeling (IPM; Chandler & Clark, 2014;Schaub & Abadi, 2011;Wilson et al., 2016), which uses multiple data sources in a joint likelihood to estimate abundance and additional demographic parameters. IPMs efficiently use data, provide a means of estimating uncertainty that is propagated among model parameters, and have the capacity to infer latent parameters for which data are not available. IPMs are versatile in the types of information they can incorporate, including capture-mark-recapture, point count, productivity, dead recovery, and telemetry data (Bled et al., 2017;Fay et al., 2019;Freeman & Crick, 2003;Horne et al., 2019;Lee et al., 2015), with the key requirement that one or more parameters are shared among the components of the IPM (Zipkin & Saunders, 2018).
Lincoln estimators (Lincoln, 1930) are increasingly used to estimate abundance at large scales (Alisauskas et al., 2009) by combining information on the rate and total number of individuals harvested; data that are typically obtained from band recoveries (Roberts et al., 2021) and hunter survey or harvest reporting. While historically underutilized, Lincoln estimators have been applied with great success in the management of multiple game species Hagen et al., 2014;Otis, 2006;Saunders et al., 2019), most notably for the harvest of waterfowl across North America (Alisauskas et al., 2014;Arnold et al., 2018). In an IPM framework, temporal dynamics in abundance can be represented in a statespace approach, and additional data sources beyond those generated by harvest may be included to inform demographic parameters estimated by the model when available (Hostetler & Chandler, 2015;Tavecchia et al., 2009).
Despite the advantages of IPMs, violating assumptions of component models can lead to bias, both for parameters directly estimated from data and those being inferred indirectly (Riecke et al., 2019), making it important to both identify potential violations and implement reasonable solutions. A major assumption of the Lincoln estimator, shared by all band recovery models, is that banding and recoveries are representative of the whole population (Alisauskas et al., 2009), and this assumption will be violated if major sources of heterogeneity in survival or harvest rates are not incorporated into the model (Pollock & Raveling, 1982). Within a contiguous population, harvest rates may vary spatially according to variable harvest regulations, land access, weather, or land cover characteristics (Burke et al., 2019;Hansen et al., 1986;Norton et al., 2012). Similarly, survival may be linked to spatially varying factors like habitat, predation risk, or weather (Fleskes et al., 2007;Perkins et al., 1997;Tolon et al., 2009). Assumptions of constant harvest rate or survival may therefore be violated across large spatial scales, which will bias estimates at finer scales. When estimating parameters statewide and applying them to management objectives that are region-specific, it is often unrealistic to assume that no heterogeneity exists among regions. Therefore, accounting for spatial heterogeneity is important to ensure accurate abundance estimates on which management will be based.
When incorporating spatial variation into models, accounting for multiple interacting factors can be difficult when each varies independently (Viana et al., 2013), and may be impractical to measure.
It may be simpler to ignore specific causes of the spatial relationship, and instead take advantage of underlying spatial correlations in the data to map the spatial structure of parameters being estimated (Cressie, 2015). Locations that are closer together in space are more likely to be similar than those farther apart (Burrough, 1995), which can facilitate covariance functions to describe the spatially-dynamic nature of a parameter. One such approach, spatial predictive processes (SPP; Banerjee et al., 2008), projects the underlying correlation among sampling sites onto a set of evenly spaced spatial knots distributed across an area of interest. SPP was initially intended as a dimension reduction approach to reduce computation requirements in Kriging for larger data sets (Banerjee et al., 2008), but the underlying framework has advantages beyond computational efficiency. For one, the covariance function does not require additional information beyond the locations of data, meaning that identifying and measuring explanatory covariates is unnecessary to represent underlying spatial heterogeneity in the process. Additionally, even spacing of spatial knots uniformly covers the area of interest, which is sometimes impractical for observed data when sampling depends on the presence of animals. Thus, using parameter estimates at evenly spaced knots may be more representative than those from sampling sites. Use of SPP is sparse within the ecological literature-although K E Y W O R D S band recovery, integrated population model, Lincoln estimator, spatial heterogeneity, spatial predictive process

T A X O N O M Y C L A S S I F I C A T I O N
Population ecology see examples for applying such an approach to estimating the spatial distribution of fisheries discards (Viana et al., 2013), avian communities (Jarzyna et al., 2014), or invasive plants (Latimer et al., 2009)-despite its broadly applicable approach for assessing spatial variation in vital rates, especially when causes of the spatial relationship are difficult to determine.
Here, we develop and present an SPP approach for incorporating spatial variation in harvest rates into an IPM to derive robust estimates of wild turkey (Meleagris gallapavo; hereafter turkey) abundance. We integrated band recovery, telemetry, and total harvest data to estimate the region-specific abundance of a two-aged (adult and juvenile), harvested population at the beginning of the hunting season ( Figure 1). Band recoveries were used to estimate harvest rate and survival rate under a modified Brownie parameterization of the dead recovery model (Brownie, 1985), in which we incorporated an SPP (Banerjee et al., 2008) to account for spatial variation in harvest rate among capture sites, and to allow estimation of harvest rates in areas where banding did not occur. To control for mortalities that occurred between banding and the beginning of the hunting season (Buderman et al., 2014), we linked survival in the band recovery model to a weekly survival rate estimated from telemetry data using a nest survival framework (Dinsmore et al., 2002). Final abundance estimates were generated using a Lincoln estimator within a statespace approach (Alisauskas et al., 2014;Lincoln, 1930). We simulated data sets under different configurations of spatial variation in harvest rate and sampling characteristics to assess bias in parameter estimates and applied the model to data collected from a field study of wild turkeys to estimate abundance and inform harvest decisions across the state of Maine.

| Band recovery model
We used a modified version of the Brownie parameterization for dead recoveries (Brownie, 1985), where recovery rates were estimated as the combined probability that a bird was killed, retrieved, and reported by a hunter within a given hunting season. We considered recovery rate synonymous with "harvest rate" (h), which was estimated as a proportion of all mortalities (harvest plus nonharvest) that occurred throughout the year. We use the term "harvest" to refer to the total number of individuals harvested and reported to MDIFW (H). We assumed 100% reporting of harvested birds, although incomplete reporting could be incorporated with additional information on reporting rate. One assumption of all band recovery models is that no mortalities occur between capture and the beginning of the first hunting season postbanding, which is likely to be violated as time between the two events increases (Cooch et al., 2021).
To better identify mortality that occurred between capture and the first hunting season, we separated each year within the conventional band recovery encounter history into two distinct and alternating occasion types, capture and the spring hunting season. The initial observation occurred during a capture occasion, and the terminal observation during a harvest occasion, resulting in two occasions per year. Survival was then differentiated according to three seasons within a year that corresponded with the intervals between recovery occasions; the period from the mid-point of capture to the first day of harvest, the interval where harvest occurs (i.e., the hunting season), and the period from last day of the harvest season to the mid-point of the following year's capture. For a given interval within the model, a Bernoulli random variable (ψ) was used to determine the probable latent survival state (z) of individual i at occasion t, where ψ was the probability of surviving all risks unrelated to reported harvests during the spring hunting season (S i,t ) given that the bird was alive at the end of the previous occasion (w i,t-1 = 1). We modeled the probability of observed data using a Bernoulli random variable (γ), Band recovery, telemetry, and total harvest information in combination with a spatial predictive process can be used in a Lincoln estimator to produce estimates of harvest rate, survival, and abundance. The flow of data and parameter estimates through the integrated population model is depicted here by a directed acyclic graph.
where y was the observed harvest of a banded individual and h was the probability a bird was harvested and reported, given that it survived all other mortality risks since the previous occasion (z i,t = 1).
Since harvests cannot occur during capture periods, we restricted h = 0 in capture occasions.
We modeled variation in h using a log-log-linear model, where β were coefficients that describe variation according to age of the individual (Juvenile, <1 year old, or Adult, >1 year old) and year of harvest, respectively.
To account for spatial variation in harvest rates, we included a mean-zero SPP (ω[c]; Viana et al., 2013; see section 2.1.2) that depended on the capture location for all individuals, c = {c 1 , c 2 , …, c n }.
Remaining variation was modeled by a nonspatial error term (ε i ), where 2.1.2 | Spatial predictive process Following the methods of Viana et al. (2013), we accounted for spatial variation in harvest rate by incorporating the mean-zero SPP into the logit-linear regression. We specified a set of evenly distributed spatial knots, c* = {c* 1 , c* 2 , …, c* m }, across the study area on which we defined a Gaussian process with exponential covariance, where σ s 2 was the spatial random effect variance and ρ was an autocorrelation function with where |d a,b | was the distance between locations c a and c b , and ϕ determined the rate of decay in correlation as distance increased between locations. To project the Gaussian process from the knots back onto harvest locations, we used the correction for bias in ε i proposed by Finley et al. (2009). Given a generic covariance function between two We then applied the covariance functions to the supplied sets of capture sites and spatial knots, which yielded site-specific estimates of the harvest rate for each capture location.

| Weekly survival rate
Weekly survival rates (s) were estimated under a nest survival modeling framework (Dinsmore et al., 2002), in which we modeled whether an individual was observed alive since its previous telemetry observation (x) as a Bernoulli random variable (μ), where s was the probability of surviving 1 week and k was the number of weeks since an individual was last observed alive. We modeled variation in s using a log-log linear model (Ergon et al., 2018), where α represented individual, temporal, and spatial regression co-

| State-space abundance estimation
Abundance (N) was derived using a Lincoln estimator (Lincoln, 1930;Alisauskas et al., 2014) for each region and timestep as where r was the region of the study area, and t was the year for which abundance was estimated. We linked abundance through time using a state-space approach. H r,t for each region were drawn from a binomial distribution where ĥ r was the mean harvest rate across years for region r and was estimated as the average harvest rate at all spatial knots (h* r ) within a region's boundaries, such that We assumed each region was closed to immigration and emigration, such that total abundance at the beginning of the hunting season (N̂a ,r,t + 1 ) was equal to the number of adults that survived the previous year (M a,r,t) , combined with juveniles that survived from the previous year and graduated to adulthood (M j,r,t ), each of which was drawn from binomial distributions where Q r was the total probability of survival, estimated as Harvest is often the most significant source of mortality for male turkeys and is dramatically greater than nonharvest mortality during the hunting season (Humberg et al., 2009); thus, we assumed that 1ĥ r represented the probability of surviving the spring hunting season by virtue of not being harvested and that nonharvest mortality during this interval was negligible. As we used categorical covariates to describe regional differences in survival and therefore could not directly estimate survival in regions in which we did not band, we estimated a

| Model validation
To assess model accuracy, we simulated data that spanned a series of adjacent regions with variable abundance. We generated a 100 km × 100 km virtual study area evenly divided into 25 regions ( Figure S1). Capture sites were randomly distributed across the study area, and each was randomly assigned either a high, medium, or low number of captured individuals. For each data set, we simulated banding, telemetry, and total harvest data for a given population, and used constant intercepts and beta coefficients to simulate weekly survival and harvest rates across the area, with modifications as described below to incorporate spatial heterogeneity. To prevent unrealistic population growth, we restricted the maximum regional abundance using a fixed carrying capacity of 5000 individuals. To introduce spatial variation into harvest rate and survival parameters, we generated Gaussian random fields using the "gstat" package (Pebesma, 2004) in program R (R Core Team, 2020), which created a location-specific beta coefficient that described spatial variation across the study area. We assessed the accuracy of estimates under variable sampling within a region and across the study area. To ensure that the simulation accurately presented a range of possible spatial heterogeneities and that the model was robust to those ranges of variation, we simulated multiple spatial configurations of harvest rate using low, medium, or high values for the partial sill, range, and nugget of the variogram used to generate the random field. In practice, this allowed us to vary the magnitude of variation in harvest rate, the maximum distance of autocorrelation, and the amount of small-scale variation in harvest rates, respectively ( Figure S2). To evaluate how the inclusion of an SPP to account for spatial heterogeneity in harvest rates impacted model estimates, we repeated the above analysis using a second model, which assumed a constant harvest rate across the study area (i.e., did not include an SPP). We then compared the relative bias between the two models.

| Case study: Wild turkeys in Maine
To demonstrate the applicability of the model, we used data collected from wild turkeys in Maine, USA. Maine is a large state (~91,647 km 2 ) with a variety of intermixed land use types and variable hunter densities. As such, we expected significant heterogeneity in harvest rates of turkeys within the state, as has been observed for other states' turkey populations (Stevens et al., 2020). Turkeys were captured during the winters (December through March) of 2018 through 2020 using rocket and drop nets and aged as either adult (>1 year old) or juvenile (<1 year old) according to plumage (Dickson, 1992). We marked turkeys with at least 1 of 4 different marking methods with associated identification numbers, including aluminum butt end leg bands, aluminum rivet bands, plastic colored leg bands, or patagial wing tags.
Nearly all individuals received at least 2 marks, and we assume retention of at least 1 mark was 100%. In addition to identification numbers, leg bands included contact information (toll-free phone number Instead, a survey was conducted to gauge success and estimate the total turkey harvest. We adjusted the model to treat 2020 total harvest data as a random variable with initial values equal to the estimated harvest from survey data. We do note that Maine also has a fall either-sex hunting season, but harvest rates were generally much lower than spring, and we received an insufficient number of recoveries to build this into the model explicitly. Therefore survival reflects all mortality (nonharvest + fall harvest) occurring between the end of the spring hunting season in year t and the beginning of the subsequent spring season in t + 1.
In Maine, turkeys are managed within 29 discrete WMDs. For the distribution of spatial knots in the SPP, we used a grid with 24 km spacing, with additional knots placed at the geographic center of each WMD. We eliminated knots from WMDs, where turkey densities (<100) and harvests (<10) were expected to be insufficient to produce reliable estimates. These knot specifications ensured that each WMD of interest had at least one knot within its boundaries and that we did not predict harvest rates beyond where our data could reasonably be considered representative.
The distance of 24 km was chosen by running multiple iterations of the model using various grid spacing. We then compared WMDspecific harvest rate estimates for each iteration to a model with a categorical covariate for WMD. We used root mean squared error to select the largest grid spacing that minimized error while also reaching convergence within the model. We compared estimates of harvest rate with and without the inclusion of the SPP component to assess whether failing to account for spatial variation in harvest rates affected parameter estimates.
Regression coefficients were given vague uniform priors. Simulation models were allowed 10,000 iterations, discarding the first 5000 from calculations. The model fit to wild turkey data was allowed 50,000 iterations, discarding the first 20,000. juveniles. We did not observe any relationship between relative bias in abundance estimates and configuration of spatial variation in harvest rates ( Figure S3). Similarly, we did not observe any differences in relative bias in abundance associated with sampling intensity for the sample sizes we considered, both for sampling within a region and for sampling across a study area within a simulation ( Figure S4).

| Simulation accuracy
When we compared relative bias in abundance as it related to the portion of the population that was banded, we found that bias became more negative as the proportion banded increased and variance increased as the proportion banded decreased ( Figure S5).
When the SPP was omitted from the model, the average relative bias for abundance estimates was 0.16 (SD = 0.14) for adults and 0.21 (SD = 0.13) for juveniles (Figure 2c,d, Figure S6)

| Case study
We captured and marked 408 male wild turkeys (187 adults We found that the estimated abundance of adult turkeys was consistently lower and had narrower credible intervals, when the SPP was F I G U R E 2 Relative bias for estimates of adult (a,c) and juvenile (b,d) abundance produced by a Lincoln estimator was on average less for a model with a spatial predictive process to account for spatial variation in harvest rates (a,b) compared with a model without (c,d). Number of simulations are depicted according to relative bias in region-specific estimates of abundance and real abundance within the region. Results shown are grouped across simulated data sets, which varied in the underlying spatial heterogeneity in the harvest rate. Estimates of adult harvest rates were consistently higher when the model included an SPP compared with the model that assumed a constant harvest rate (Figure 6c). For juveniles, we observed substantially less difference in parameter estimates between models with and without an SPP (Figure 6b,d), consistent with a lower overall harvest rate for juvenile males with inherently less room for variability as a result.

| DISCUSS ION
We found that a Lincoln estimator incorporating an SPP in the band recovery model consistently estimated abundance with low bias across a variety of spatial configurations of harvest rate for simulated data. When the SPP was omitted from the model, relative bias in abundance increased for all configurations of spatial heterogeneity in harvest rate. Additionally, we found no difference in relative bias according to either sampling intensity or the underlying nature of spatial heterogeneity in harvest rate. When applied to real data collected from turkeys in Maine, we observed a wide range of harvest rates and abundance among wildlife management districts. As expected, the basic model that did not incorporate spatial heterogeneity in harvest rates via the SPP underestimated harvest rates in some WMDs, which resulted in an overestimation of abundance in those districts, and statewide. The variation in harvest rates we observed is typical for wild turkey populations and should be expected especially across large spatial scales (Norton et al., 2012). Most recent applications of Lincoln estimators have treated harvest rates as uniform across large areas (Alisauskas et al., 2014;Hagen et al., 2018;Shirkey & Gates, 2020). While this assumption may sometimes be appropriate and will likely depend on the definition of "large areas," our results suggest that future applications could consider more explicitly incorporating spatial variation in harvest rates to improve inference. While we chose to use an SPP to accomplish this, other methods such as simultaneous autoregressive and conditional autoregressive models may also be a viable option (Fortin et al., 2012).
Harvest management decisions must often consider populations that span considerably large spatial scales (Robinson et al., 2016), and it may not be necessary to reconcile inference at the local scale when decisions are made at regional levels (Johnson et al., 2015). Instead, estimates can be aggregated to summarize relationships within a region's boundaries, making fine-scale differences in parameters at local scales less important than adequately capturing the general trend in a parameter across space. Aggregating estimates to describe parameters by region can be performed using multiple methods, with the simplest solution being to average estimates within each region.
However, consideration must be given to the sampling design used, as clumped or sparse sampling within an area could lead to bias if sampled locations differ greatly from the mean across a landscape (Hooten et al., 2017). To some degree, SPP can mitigate such issues by using the entire data set to define a spatial correlation function and projecting it onto evenly distributed spatial knots from which estimates are then made. This will have the benefit of smoothing the prediction surface, minimizing the impact of any single sampling location, which may otherwise have outsized impacts on local averages. However, if larger areas or districts have particularly high or low harvest rates, and are unsampled, the model would not be able to interpolate those relationships, meaning that adequate and representative sampling is still an important component of study design.
When adequate sampling is performed, SPP has proven to be effective at identifying localized "hotspots" in variation (Viana et al., 2013), and indeed we found that model predictions were robust to a wide range of underlying spatial heterogeneity in harvest rate.
We observed some variation in the magnitude and direction of error across iterations, which is common when assessing IPMs (Abadi et al., 2010;Fieberg et al., 2010). We further found that relative bias became more negative as the proportion of a population that was banded increased. This shift in bias is consistent with a long-understood relationship, where the accuracy of markrecapture models depends on the ratio of the banded sample to total abundance (Robson & Regier, 1964), and reinforces previous recommendations that sample size of banding studies should be informed by expected population size (Robson & Regier, 1964).

F I G U R E 3
Spatial knots (red triangles) chosen for the spatial predictive process provide more uniform and complete coverage of the area of interest compared with capture sites (blue dots). Spatial knots were distributed across a 24 km × 24 km grid with additional knots placed at the geographic center of each WMD. We eliminated knots from WMDs where turkey densities (<100) and harvests (<10) were expected to be insufficient to produce reliable estimates. Sample size for capture sites is indicated by the size of the dot, with large dots meaning larger sample sizes. Color of dots indicates whether a site had telemetry devices deployed (dark blue) or did not (light blue).
Despite the advantages of SPP, there are still opportunities for improvement. Although SPP should be notably faster than alternatives, especially as the amount of data increases, using MCMC can still lead to lengthy computing times (Banerjee et al., 2008) so alternative posterior samplers should be considered to decrease processing time.
The advantages of SPP will not overcome extremely sparse data availability or poor sampling design. As previously mentioned, the number of individuals marked should be proportional to the expected population size being sampled. While we did not observe an effect of sample size on the model bias for the sample sizes we considered, further exploration with more limited data sets may be necessary to find a threshold at which estimates are no longer reliable. Similarly, the distribution of sampling locations and the configuration of spatial knots across a study area should be informed by the ecology of the system being studied. The number and placement of knots for the SPP are not trivial, and while there appears to be a wide margin for error, these decisions have an impact on estimates (Banerjee et al., 2008).

F I G U R E 5
Harvest rate estimates were variable among wildlife management districts across Maine, USA, for both male and juvenile wild turkeys. Regional differences in harvest rates are depicted for adult male turkeys with color indicating the mean harvest rate for each management district.
experiencing spatial heterogeneity in one parameter could result in increased variation in estimates of the other.
Multiple methods are currently used to monitor turkey populations. Many states approximate population trends using spring harvest data (e.g., CDEEP, 2016; Harms et al., 2017;Healy, 2000), but this does not provide an estimate of true abundance, which may be preferable for setting regulations (Lint et al., 1995). This method also requires accounting for changes in the abundance of birds and harvest rates, both of which influence the number of birds harvested through time (Paloheimo & Fraser, 1981). Surveys such as summer sighting (PGC, 2021), gobble counts (Rioux et al., 2009), and camera traps (Gonnerman, 2017) can be used to produce estimates of population size at smaller scales but are unrealistic to implement for statewide management. An IPM, such as the one we have implemented, provides a data-driven alternative that can be scaled to the scope of turkey management decisions. It is relatively costeffective as it uses often already implemented mandatory reporting of harvests and only requires periodic captures of individuals for F I G U R E 6 Estimates of adult abundance and harvest rate differed between Lincoln estimators with and without a spatial predictive process included, but juvenile harvest rate and abundance were largely similar. Estimates from both models are presented with associated error bars for abundance (a,b; shown on the log scale) and harvest rate (c,d). For all figures, circles correspond to WMD-by-year estimates. For c and d, triangles indicate the harvest rate estimated without the inclusion of the SPP, and circles indicate the SPP.
banding. Similar IPMs have been implemented for turkey and waterfowl populations to great success (Arnold et al., 2018;Diefenbach et al., 2012), demonstrating that this is a feasible alternative that, with the inclusion of an SPP component, overcomes many of the shortcomings of more common monitoring methods. We imagine this approach may have more broad applicability to other similarlymanaged harvest systems, such as those for large mammals, where regional or subregional variation in harvest regulations is common.

| CON CLUS IONS
Management decisions based on biased estimates of abundance may lead to harvest regulations that exceed sustainable levels or are unnecessarily restrictive (Dillingham & Fletcher, 2008). Violations of the Lincoln estimator's assumption of representative harvest may result in such bias. As Lincoln estimators become more widely applied, it is important to consider a mechanism to account for spatial variation in harvest rates, which can vary according to a broad range of spatially varying ecological, environmental, and socio-economic factors that are sometimes difficult to measure (Pope & Powell, 2021). For such cases, the combination of Lincoln estimator and SPP is an especially relevant tool for capturing the magnitude and distribution of variation to reduce bias in estimates used for management. SPP functions as a component of a generalized linear mixed model framework, making it compatible with many analytical methods currently used in ecology, and therefore should be accessible to those who are less familiar with spatial statistics. While we chose to apply these methods to harvest rates within a band recovery model, the use of SPP should be widely applicable across methods for vital rate estimation. writing -original draft (equal); writing -review and editing (lead).

ACK N OWLED G EM ENTS
We thank R. B. Allen, all the Regional Wildlife Biologists, many other MDIFW staff, and especially B. Currier, B. Peterson, K. Leary, G.
LeClaire, K. Overturf, and N. Aielo, who each made significant contributions to field data collection. Student volunteers from the University of Maine and Unity College also contributed to field work. We are grateful to Dr. P. Milligan for providing lab protocols, primers, and positive controls and B. Tweedie and C. Desjardins for their contributions in the lab.
We thank American Forest Management and the numerous individual private owners for land access and accommodations. This is Maine Agricultural and Forest Experiment Station Publication no. 3791.

CO N FLI C T O F I NTE R E S T
All authors involved attest to having no conflicts of interest related to the submitted manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in GitHub at https://github.com/mattg onner man/BandR ecove ryModel (Gonnerman, 2022).